home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dgebrd.f < prev    next >
Text File  |  1996-07-19  |  9KB  |  259 lines
  1.       SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
  2.      $                   INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       INTEGER            INFO, LDA, LWORK, M, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   A( LDA, * ), D( * ), E( * ), TAUP( * ),
  14.      $                   TAUQ( * ), WORK( LWORK )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  DGEBRD reduces a general real M-by-N matrix A to upper or lower
  21. *  bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
  22. *
  23. *  If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
  24. *
  25. *  Arguments
  26. *  =========
  27. *
  28. *  M       (input) INTEGER
  29. *          The number of rows in the matrix A.  M >= 0.
  30. *
  31. *  N       (input) INTEGER
  32. *          The number of columns in the matrix A.  N >= 0.
  33. *
  34. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  35. *          On entry, the M-by-N general matrix to be reduced.
  36. *          On exit,
  37. *          if m >= n, the diagonal and the first superdiagonal are
  38. *            overwritten with the upper bidiagonal matrix B; the
  39. *            elements below the diagonal, with the array TAUQ, represent
  40. *            the orthogonal matrix Q as a product of elementary
  41. *            reflectors, and the elements above the first superdiagonal,
  42. *            with the array TAUP, represent the orthogonal matrix P as
  43. *            a product of elementary reflectors;
  44. *          if m < n, the diagonal and the first subdiagonal are
  45. *            overwritten with the lower bidiagonal matrix B; the
  46. *            elements below the first subdiagonal, with the array TAUQ,
  47. *            represent the orthogonal matrix Q as a product of
  48. *            elementary reflectors, and the elements above the diagonal,
  49. *            with the array TAUP, represent the orthogonal matrix P as
  50. *            a product of elementary reflectors.
  51. *          See Further Details.
  52. *
  53. *  LDA     (input) INTEGER
  54. *          The leading dimension of the array A.  LDA >= max(1,M).
  55. *
  56. *  D       (output) DOUBLE PRECISION array, dimension (min(M,N))
  57. *          The diagonal elements of the bidiagonal matrix B:
  58. *          D(i) = A(i,i).
  59. *
  60. *  E       (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
  61. *          The off-diagonal elements of the bidiagonal matrix B:
  62. *          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
  63. *          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
  64. *
  65. *  TAUQ    (output) DOUBLE PRECISION array dimension (min(M,N))
  66. *          The scalar factors of the elementary reflectors which
  67. *          represent the orthogonal matrix Q. See Further Details.
  68. *
  69. *  TAUP    (output) DOUBLE PRECISION array, dimension (min(M,N))
  70. *          The scalar factors of the elementary reflectors which
  71. *          represent the orthogonal matrix P. See Further Details.
  72. *
  73. *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  74. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  75. *
  76. *  LWORK   (input) INTEGER
  77. *          The length of the array WORK.  LWORK >= max(1,M,N).
  78. *          For optimum performance LWORK >= (M+N)*NB, where NB
  79. *          is the optimal blocksize.
  80. *
  81. *  INFO    (output) INTEGER
  82. *          = 0:  successful exit
  83. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  84. *
  85. *  Further Details
  86. *  ===============
  87. *
  88. *  The matrices Q and P are represented as products of elementary
  89. *  reflectors:
  90. *
  91. *  If m >= n,
  92. *
  93. *     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)
  94. *
  95. *  Each H(i) and G(i) has the form:
  96. *
  97. *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
  98. *
  99. *  where tauq and taup are real scalars, and v and u are real vectors;
  100. *  v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i);
  101. *  u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n);
  102. *  tauq is stored in TAUQ(i) and taup in TAUP(i).
  103. *
  104. *  If m < n,
  105. *
  106. *     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)
  107. *
  108. *  Each H(i) and G(i) has the form:
  109. *
  110. *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
  111. *
  112. *  where tauq and taup are real scalars, and v and u are real vectors;
  113. *  v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
  114. *  u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
  115. *  tauq is stored in TAUQ(i) and taup in TAUP(i).
  116. *
  117. *  The contents of A on exit are illustrated by the following examples:
  118. *
  119. *  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
  120. *
  121. *    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
  122. *    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
  123. *    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
  124. *    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
  125. *    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
  126. *    (  v1  v2  v3  v4  v5 )
  127. *
  128. *  where d and e denote diagonal and off-diagonal elements of B, vi
  129. *  denotes an element of the vector defining H(i), and ui an element of
  130. *  the vector defining G(i).
  131. *
  132. *  =====================================================================
  133. *
  134. *     .. Parameters ..
  135.       DOUBLE PRECISION   ONE
  136.       PARAMETER          ( ONE = 1.0D+0 )
  137. *     ..
  138. *     .. Local Scalars ..
  139.       INTEGER            I, IINFO, J, LDWRKX, LDWRKY, MINMN, NB, NBMIN,
  140.      $                   NX
  141.       DOUBLE PRECISION   WS
  142. *     ..
  143. *     .. External Subroutines ..
  144.       EXTERNAL           DGEBD2, DGEMM, DLABRD, XERBLA
  145. *     ..
  146. *     .. Intrinsic Functions ..
  147.       INTRINSIC          MAX, MIN
  148. *     ..
  149. *     .. External Functions ..
  150.       INTEGER            ILAENV
  151.       EXTERNAL           ILAENV
  152. *     ..
  153. *     .. Executable Statements ..
  154. *
  155. *     Test the input parameters
  156. *
  157.       INFO = 0
  158.       IF( M.LT.0 ) THEN
  159.          INFO = -1
  160.       ELSE IF( N.LT.0 ) THEN
  161.          INFO = -2
  162.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  163.          INFO = -4
  164.       ELSE IF( LWORK.LT.MAX( 1, M, N ) ) THEN
  165.          INFO = -10
  166.       END IF
  167.       IF( INFO.LT.0 ) THEN
  168.          CALL XERBLA( 'DGEBRD', -INFO )
  169.          RETURN
  170.       END IF
  171. *
  172. *     Quick return if possible
  173. *
  174.       MINMN = MIN( M, N )
  175.       IF( MINMN.EQ.0 ) THEN
  176.          WORK( 1 ) = 1
  177.          RETURN
  178.       END IF
  179. *
  180.       WS = MAX( M, N )
  181.       LDWRKX = M
  182.       LDWRKY = N
  183. *
  184. *     Set the block size NB and the crossover point NX.
  185. *
  186.       NB = MAX( 1, ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) )
  187. *
  188.       IF( NB.GT.1 .AND. NB.LT.MINMN ) THEN
  189. *
  190. *        Determine when to switch from blocked to unblocked code.
  191. *
  192.          NX = MAX( NB, ILAENV( 3, 'DGEBRD', ' ', M, N, -1, -1 ) )
  193.          IF( NX.LT.MINMN ) THEN
  194.             WS = ( M+N )*NB
  195.             IF( LWORK.LT.WS ) THEN
  196. *
  197. *              Not enough work space for the optimal NB, consider using
  198. *              a smaller block size.
  199. *
  200.                NBMIN = ILAENV( 2, 'DGEBRD', ' ', M, N, -1, -1 )
  201.                IF( LWORK.GE.( M+N )*NBMIN ) THEN
  202.                   NB = LWORK / ( M+N )
  203.                ELSE
  204.                   NB = 1
  205.                   NX = MINMN
  206.                END IF
  207.             END IF
  208.          END IF
  209.       ELSE
  210.          NX = MINMN
  211.       END IF
  212. *
  213.       DO 30 I = 1, MINMN - NX, NB
  214. *
  215. *        Reduce rows and columns i:i+nb-1 to bidiagonal form and return
  216. *        the matrices X and Y which are needed to update the unreduced
  217. *        part of the matrix
  218. *
  219.          CALL DLABRD( M-I+1, N-I+1, NB, A( I, I ), LDA, D( I ), E( I ),
  220.      $                TAUQ( I ), TAUP( I ), WORK, LDWRKX,
  221.      $                WORK( LDWRKX*NB+1 ), LDWRKY )
  222. *
  223. *        Update the trailing submatrix A(i+nb:m,i+nb:n), using an update
  224. *        of the form  A := A - V*Y' - X*U'
  225. *
  226.          CALL DGEMM( 'No transpose', 'Transpose', M-I-NB+1, N-I-NB+1,
  227.      $               NB, -ONE, A( I+NB, I ), LDA,
  228.      $               WORK( LDWRKX*NB+NB+1 ), LDWRKY, ONE,
  229.      $               A( I+NB, I+NB ), LDA )
  230.          CALL DGEMM( 'No transpose', 'No transpose', M-I-NB+1, N-I-NB+1,
  231.      $               NB, -ONE, WORK( NB+1 ), LDWRKX, A( I, I+NB ), LDA,
  232.      $               ONE, A( I+NB, I+NB ), LDA )
  233. *
  234. *        Copy diagonal and off-diagonal elements of B back into A
  235. *
  236.          IF( M.GE.N ) THEN
  237.             DO 10 J = I, I + NB - 1
  238.                A( J, J ) = D( J )
  239.                A( J, J+1 ) = E( J )
  240.    10       CONTINUE
  241.          ELSE
  242.             DO 20 J = I, I + NB - 1
  243.                A( J, J ) = D( J )
  244.                A( J+1, J ) = E( J )
  245.    20       CONTINUE
  246.          END IF
  247.    30 CONTINUE
  248. *
  249. *     Use unblocked code to reduce the remainder of the matrix
  250. *
  251.       CALL DGEBD2( M-I+1, N-I+1, A( I, I ), LDA, D( I ), E( I ),
  252.      $             TAUQ( I ), TAUP( I ), WORK, IINFO )
  253.       WORK( 1 ) = WS
  254.       RETURN
  255. *
  256. *     End of DGEBRD
  257. *
  258.       END
  259.